chi22-course
08 March, 2022
Introduction
In this document, we will outline the Bayesian analogs of the statistical analyses described here (and in the previous course lecture).
Helper functions.
# this function extracts results from different models and generate results of the same format to be used in visualizations
tidy.wrapper = function(model) {
if (class(model) == "lm") {
tidy(model, conf.int = TRUE) %>%
select(-c(statistic, p.value)) %>%
mutate(model = "Frequentist") %>%
select(model, everything())
} else if (class(model) == "brmsfit") {
tidy(model) %>%
filter(effect == "fixed") %>%
select(c(term:conf.high)) %>%
mutate(model = "Bayesian") %>%
select(model, everything())
} else {
stop("unknown model class")
}
}TODO: explain the dataset here
Dataset 1
Let’s first load and look at the data in a table
# Fumemng -> Abhraneel I'm using the Rproj in the root folder, although my path is identical...
dataset = readr::read_csv("02-bayesian_stats/data/blinded.csv")##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## experiment = col_double(),
## condition = col_character(),
## effectiveness = col_double()
## )# dataset = readr::read_csv("~/Documents/Github/chi22-course/data/blinded.csv")
kable(head(dataset))| experiment | condition | effectiveness |
|---|---|---|
| 1 | no_graph | 7 |
| 1 | graph | 6 |
| 1 | no_graph | 9 |
| 1 | no_graph | 4 |
| 1 | graph | 6 |
| 1 | no_graph | 7 |
dataset %>%
mutate(effectiveness = fct_rev(factor(effectiveness, levels = 1:9))) %>%
# stacked bar plot
ggplot(aes(x = condition, fill = effectiveness)) +
geom_bar(position = "stack", stat="count") +
# plot data for different experiments as small multiples
facet_wrap( ~ experiment) +
# grey color scale is robust for colorblind
scale_fill_brewer(palette="Greys", drop = FALSE) +
# horizontal plot
coord_flip() +
# legend
guides(fill = guide_legend(reverse = TRUE)) 
As we can see above, the original dataset contains results from four different experiments. For the purposes of this lecture, we will confine ourselves to the first experiment.
exp1.data = dataset %>%
filter(experiment == 1)
exp1.data## # A tibble: 123 × 3
## experiment condition effectiveness
## <dbl> <chr> <dbl>
## 1 1 no_graph 7
## 2 1 graph 6
## 3 1 no_graph 9
## 4 1 no_graph 4
## 5 1 graph 6
## 6 1 no_graph 7
## 7 1 graph 7
## 8 1 no_graph 7
## 9 1 no_graph 7
## 10 1 graph 6
## # … with 113 more rowsModel 1. Wilcoxon signed rank test
This is a non-parametric test which we will skip for now. Although, there exists Bayesian non-parametric methods, they are more advanced for this lecture.
Model 2. t-test
This is the linear model equivalent for the paired sample t-test
dataset1.lm.freqt <-
lm(
effectiveness ~ condition - 1,
data = exp1.data
) exp1.data %>%
ggplot(., aes(x = effectiveness)) +
geom_histogram(fill = '#351c75', color = NA)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
dataset1.brm.bayesiant.formula <- bf(
effectiveness ~ condition,
family = student()
)kable(as.tibble(get_prior(dataset1.brm.bayesiant.formula, data = exp1.data)))## Warning: `as.tibble()` was deprecated in tibble 2.0.0.
## Please use `as_tibble()` instead.
## The signature and semantics have changed, see `?as_tibble`.| prior | class | coef | group | resp | dpar | nlpar | bound | source |
|---|---|---|---|---|---|---|---|---|
| b | default | |||||||
| b | conditionno_graph | default | ||||||
| student_t(3, 7, 2.5) | Intercept | default | ||||||
| gamma(2, 0.1) | nu | default | ||||||
| student_t(3, 0, 2.5) | sigma | default |
dataset1.brm.bayesiant.priors = c(
prior(normal(0, 1), class = "b"), # there's a lot of data so even fairly "strong" priors are going to not matter so much here
prior(student_t(3, 0, 1), class = "sigma")
)dataset1.brm.bayesiant.priorchecks <-
brm(
dataset1.brm.bayesiant.formula,
prior = dataset1.brm.bayesiant.priors,
data = exp1.data,
backend = "cmdstanr",
sample_prior = "only",
file = "02-bayesian_stats/rds/dataset1.brm.bayesiant.priorchecks.rds"
)## Compiling Stan program...## Start sampling## Running MCMC with 4 sequential chains...
##
## Chain 1 Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1 Iteration: 100 / 2000 [ 5%] (Warmup)
## Chain 1 Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1 Iteration: 300 / 2000 [ 15%] (Warmup)
## Chain 1 Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1 Iteration: 500 / 2000 [ 25%] (Warmup)
## Chain 1 Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1 Iteration: 700 / 2000 [ 35%] (Warmup)
## Chain 1 Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1 Iteration: 900 / 2000 [ 45%] (Warmup)
## Chain 1 Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1 Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1 Iteration: 1100 / 2000 [ 55%] (Sampling)
## Chain 1 Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1 Iteration: 1300 / 2000 [ 65%] (Sampling)
## Chain 1 Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1 Iteration: 1500 / 2000 [ 75%] (Sampling)
## Chain 1 Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1 Iteration: 1700 / 2000 [ 85%] (Sampling)
## Chain 1 Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1 Iteration: 1900 / 2000 [ 95%] (Sampling)
## Chain 1 Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2 Iteration: 100 / 2000 [ 5%] (Warmup)
## Chain 2 Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2 Iteration: 300 / 2000 [ 15%] (Warmup)
## Chain 2 Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2 Iteration: 500 / 2000 [ 25%] (Warmup)
## Chain 2 Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2 Iteration: 700 / 2000 [ 35%] (Warmup)
## Chain 2 Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2 Iteration: 900 / 2000 [ 45%] (Warmup)
## Chain 2 Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2 Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2 Iteration: 1100 / 2000 [ 55%] (Sampling)
## Chain 2 Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2 Iteration: 1300 / 2000 [ 65%] (Sampling)
## Chain 2 Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2 Iteration: 1500 / 2000 [ 75%] (Sampling)
## Chain 2 Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2 Iteration: 1700 / 2000 [ 85%] (Sampling)
## Chain 2 Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2 Iteration: 1900 / 2000 [ 95%] (Sampling)
## Chain 2 Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2 finished in 0.0 seconds.
## Chain 3 Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3 Iteration: 100 / 2000 [ 5%] (Warmup)
## Chain 3 Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3 Iteration: 300 / 2000 [ 15%] (Warmup)
## Chain 3 Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3 Iteration: 500 / 2000 [ 25%] (Warmup)
## Chain 3 Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3 Iteration: 700 / 2000 [ 35%] (Warmup)
## Chain 3 Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3 Iteration: 900 / 2000 [ 45%] (Warmup)
## Chain 3 Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3 Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3 Iteration: 1100 / 2000 [ 55%] (Sampling)
## Chain 3 Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3 Iteration: 1300 / 2000 [ 65%] (Sampling)
## Chain 3 Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3 Iteration: 1500 / 2000 [ 75%] (Sampling)
## Chain 3 Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3 Iteration: 1700 / 2000 [ 85%] (Sampling)
## Chain 3 Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3 Iteration: 1900 / 2000 [ 95%] (Sampling)
## Chain 3 Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3 finished in 0.0 seconds.
## Chain 4 Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4 Iteration: 100 / 2000 [ 5%] (Warmup)
## Chain 4 Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4 Iteration: 300 / 2000 [ 15%] (Warmup)
## Chain 4 Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4 Iteration: 500 / 2000 [ 25%] (Warmup)
## Chain 4 Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4 Iteration: 700 / 2000 [ 35%] (Warmup)
## Chain 4 Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4 Iteration: 900 / 2000 [ 45%] (Warmup)
## Chain 4 Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4 Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4 Iteration: 1100 / 2000 [ 55%] (Sampling)
## Chain 4 Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4 Iteration: 1300 / 2000 [ 65%] (Sampling)
## Chain 4 Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4 Iteration: 1500 / 2000 [ 75%] (Sampling)
## Chain 4 Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4 Iteration: 1700 / 2000 [ 85%] (Sampling)
## Chain 4 Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4 Iteration: 1900 / 2000 [ 95%] (Sampling)
## Chain 4 Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 1.0 seconds.dataset1.bayesiant.yprior <-
posterior_predict(dataset1.brm.bayesiant.priorchecks)
ggplot() +
geom_density(aes(x = dataset1.bayesiant.yprior),
color = '#351c75',
alpha = .5,
size = 1) +
xlab('prior draws') +
ggtitle('prior preditive check')
dataset1.brm.bayesiant =
brm(dataset1.brm.bayesiant.formula,
prior = dataset1.brm.bayesiant.priors,
data = exp1.data,
#You may not need this if rstan works with your computer
backend = "cmdstanr",
#Save the model
file = "02-bayesian_stats/rds/dataset1.brm.bayesiant.rds"
)## Compiling Stan program...## Start sampling## Running MCMC with 4 sequential chains...
##
## Chain 1 Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1 Iteration: 100 / 2000 [ 5%] (Warmup)
## Chain 1 Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1 Iteration: 300 / 2000 [ 15%] (Warmup)
## Chain 1 Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1 Iteration: 500 / 2000 [ 25%] (Warmup)
## Chain 1 Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1 Iteration: 700 / 2000 [ 35%] (Warmup)
## Chain 1 Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1 Iteration: 900 / 2000 [ 45%] (Warmup)
## Chain 1 Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1 Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1 Iteration: 1100 / 2000 [ 55%] (Sampling)
## Chain 1 Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1 Iteration: 1300 / 2000 [ 65%] (Sampling)
## Chain 1 Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1 Iteration: 1500 / 2000 [ 75%] (Sampling)
## Chain 1 Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1 Iteration: 1700 / 2000 [ 85%] (Sampling)
## Chain 1 Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1 Iteration: 1900 / 2000 [ 95%] (Sampling)
## Chain 1 Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2 Iteration: 100 / 2000 [ 5%] (Warmup)
## Chain 2 Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2 Iteration: 300 / 2000 [ 15%] (Warmup)
## Chain 2 Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2 Iteration: 500 / 2000 [ 25%] (Warmup)
## Chain 2 Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2 Iteration: 700 / 2000 [ 35%] (Warmup)
## Chain 2 Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2 Iteration: 900 / 2000 [ 45%] (Warmup)
## Chain 2 Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2 Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2 Iteration: 1100 / 2000 [ 55%] (Sampling)
## Chain 2 Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2 Iteration: 1300 / 2000 [ 65%] (Sampling)
## Chain 2 Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2 Iteration: 1500 / 2000 [ 75%] (Sampling)
## Chain 2 Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2 Iteration: 1700 / 2000 [ 85%] (Sampling)
## Chain 2 Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2 Iteration: 1900 / 2000 [ 95%] (Sampling)
## Chain 2 Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2 finished in 0.1 seconds.
## Chain 3 Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3 Iteration: 100 / 2000 [ 5%] (Warmup)
## Chain 3 Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3 Iteration: 300 / 2000 [ 15%] (Warmup)
## Chain 3 Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3 Iteration: 500 / 2000 [ 25%] (Warmup)
## Chain 3 Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3 Iteration: 700 / 2000 [ 35%] (Warmup)
## Chain 3 Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3 Iteration: 900 / 2000 [ 45%] (Warmup)
## Chain 3 Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3 Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3 Iteration: 1100 / 2000 [ 55%] (Sampling)
## Chain 3 Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3 Iteration: 1300 / 2000 [ 65%] (Sampling)
## Chain 3 Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3 Iteration: 1500 / 2000 [ 75%] (Sampling)
## Chain 3 Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3 Iteration: 1700 / 2000 [ 85%] (Sampling)
## Chain 3 Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3 Iteration: 1900 / 2000 [ 95%] (Sampling)
## Chain 3 Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3 finished in 0.1 seconds.
## Chain 4 Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4 Iteration: 100 / 2000 [ 5%] (Warmup)
## Chain 4 Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4 Iteration: 300 / 2000 [ 15%] (Warmup)
## Chain 4 Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4 Iteration: 500 / 2000 [ 25%] (Warmup)
## Chain 4 Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4 Iteration: 700 / 2000 [ 35%] (Warmup)
## Chain 4 Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4 Iteration: 900 / 2000 [ 45%] (Warmup)
## Chain 4 Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4 Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4 Iteration: 1100 / 2000 [ 55%] (Sampling)
## Chain 4 Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4 Iteration: 1300 / 2000 [ 65%] (Sampling)
## Chain 4 Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4 Iteration: 1500 / 2000 [ 75%] (Sampling)
## Chain 4 Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4 Iteration: 1700 / 2000 [ 85%] (Sampling)
## Chain 4 Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4 Iteration: 1900 / 2000 [ 95%] (Sampling)
## Chain 4 Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.9 seconds. summary(dataset1.brm.bayesiant)## Family: student
## Links: mu = identity; sigma = identity; nu = identity
## Formula: effectiveness ~ condition
## Data: exp1.data (Number of observations: 123)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 6.63 0.18 6.29 6.97 1.00 3579 2545
## conditionno_graph 0.16 0.24 -0.31 0.63 1.00 3592 2585
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 1.29 0.12 1.07 1.52 1.00 3264 2435
## nu 16.53 11.57 4.27 46.43 1.00 3170 2975
##
## Samples were drawn using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).color_scheme_set(scheme = "purple")
plot(dataset1.brm.bayesiant)
dataset1.bayesiant.y <- exp1.data$effectiveness
dataset1.bayesiant.yrep <- posterior_predict(dataset1.brm.bayesiant, draws = 50)
kable(head(as_tibble(dataset1.bayesiant.yrep[,1:11]))) ## Warning: The `x` argument of `as_tibble.matrix()` must have unique column names if `.name_repair` is omitted as of tibble 2.0.0.
## Using compatibility `.name_repair`.| V1 | V2 | V3 | V4 | V5 | V6 | V7 | V8 | V9 | V10 | V11 |
|---|---|---|---|---|---|---|---|---|---|---|
| 7.350071 | 7.550894 | 4.346016 | 7.402838 | 6.105783 | 6.529116 | 8.815189 | 8.225576 | 6.937635 | 7.539691 | 5.342559 |
| 7.294639 | 7.101592 | 4.382338 | 8.122888 | 9.732572 | 5.105547 | 5.986431 | 7.017659 | 7.960890 | 6.007279 | 5.592948 |
| 6.181452 | 2.912017 | 5.652890 | 6.575436 | 7.737640 | 7.520308 | 8.458674 | 6.613323 | 8.191181 | 7.433597 | 6.322437 |
| 8.535437 | 7.939905 | 5.246086 | 5.892014 | 4.370619 | 6.351047 | 6.903111 | 4.641686 | 7.422193 | 7.811355 | 5.179726 |
| 7.361354 | 4.894749 | 6.116840 | 5.985281 | 6.603010 | 6.110403 | 6.298743 | 7.107104 | 5.918865 | 6.950818 | 6.126804 |
| 6.634058 | 8.767211 | 7.708752 | 4.349791 | 6.250896 | 5.995275 | 6.298566 | 9.387987 | 5.926242 | 6.926282 | 8.637597 |
ppc_hist(y = dataset1.bayesiant.y,
yrep = dataset1.bayesiant.yrep[1:8, ],
binwidth = .5)
ppc_dens_overlay(y = dataset1.bayesiant.y,
yrep = dataset1.bayesiant.yrep[1:30, ])
ppc_stat_grouped(y = dataset1.bayesiant.y,
yrep = dataset1.bayesiant.yrep,
group = exp1.data$condition)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Compared to the frequentist estimates:
bind_rows(tidy.wrapper(dataset1.lm.freqt),
tidy.wrapper(dataset1.brm.bayesiant)) %>%
ggplot() +
geom_pointrange(
aes(
x = model,
y = estimate,
ymin = conf.low,
ymax = conf.high,
color = term
),
position = position_dodge(width = 0.2)
) +
scale_color_brewer(palette = "Set1") +
ylab('effectiveness') +
scale_y_continuous(breaks = 1:9, limits = c(1, 9)) +
coord_flip()
dataset1.bayesiant.posterior_epred <- tibble(condition = c('graph', 'no_graph')) %>%
add_epred_draws(dataset1.brm.bayesiant, re_formula = NA, allow_new_levels = FALSE) %>%
ungroup()
kable(head(dataset1.bayesiant.posterior_epred))| condition | .row | .chain | .iteration | .draw | .epred |
|---|---|---|---|---|---|
| graph | 1 | NA | NA | 1 | 6.50377 |
| graph | 1 | NA | NA | 2 | 6.64234 |
| graph | 1 | NA | NA | 3 | 6.59090 |
| graph | 1 | NA | NA | 4 | 6.50090 |
| graph | 1 | NA | NA | 5 | 6.69152 |
| graph | 1 | NA | NA | 6 | 6.47357 |
dataset1.bayesiant.posterior_comparison <- dataset1.bayesiant.posterior_epred %>%
select(-c(.chain, .iteration, .row)) %>%
compare_levels(variable = .epred, by = condition)
kable(head(dataset1.bayesiant.posterior_comparison))| .draw | condition | .epred |
|---|---|---|
| 1 | no_graph - graph | 0.1420210 |
| 2 | no_graph - graph | 0.1585310 |
| 3 | no_graph - graph | 0.3108710 |
| 4 | no_graph - graph | 0.0574165 |
| 5 | no_graph - graph | 0.3252280 |
| 6 | no_graph - graph | 0.2323100 |
kable(dataset1.bayesiant.posterior_comparison %>%
mean_qi(.epred))| condition | .epred | .lower | .upper | .width | .point | .interval |
|---|---|---|---|---|---|---|
| no_graph - graph | 0.1609358 | -0.3058675 | 0.6307221 | 0.95 | mean | qi |
dataset1.bayesiant.posterior_comparison %>%
median_qi(.epred) %>%
ggplot() +
geom_point(aes(x = .epred, y = condition), color = DARK_PURPLE, size = 3) +
geom_errorbarh(
aes(xmin = .lower, xmax = .upper, y = condition),
color = DARK_PURPLE,
alpha = .5,
size = 2,
height = 0
) +
geom_vline(aes(xintercept = 0), linetype = "dashed", color = "gray") +
coord_cartesian(ylim = c(0, 2), xlim = c(-1, 1)) +
xlab('') + ylab('')
Model 3. Ordinal Logistic Regression
dataset1.brm.olr.formula <-
bf(
effectiveness ~ condition,
family = cumulative("logit")
)get_prior(dataset1.brm.olr.formula, data = exp1.data)## prior class coef group resp dpar nlpar bound
## (flat) b
## (flat) b conditionno_graph
## student_t(3, 0, 2.5) Intercept
## student_t(3, 0, 2.5) Intercept 1
## student_t(3, 0, 2.5) Intercept 2
## student_t(3, 0, 2.5) Intercept 3
## student_t(3, 0, 2.5) Intercept 4
## student_t(3, 0, 2.5) Intercept 5
## student_t(3, 0, 2.5) Intercept 6
## student_t(3, 0, 2.5) Intercept 7
## student_t(3, 0, 2.5) Intercept 8
## source
## default
## (vectorized)
## default
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)dataset1.brm.olr.priors = c(
prior(normal(0, 1), class = "b"),
prior(student_t(3, 0, 1), class = "Intercept")
)dataset1.brm.olr.priorchecks <- brm(
effectiveness ~ condition,
prior = dataset1.brm.olr.priors,
data = exp1.data,
family= cumulative("logit"),
backend = 'cmdstanr',
sample_prior = 'only',
file = "02-bayesian_stats/rds/dataset1.brm.bayesiant.priorchecks.rds"
)
dataset1.olr.yprior <-
posterior_predict(dataset1.brm.olr.priorchecks)
ggplot() +
geom_histogram(aes(x = dataset1.olr.yprior),
fill = '#351c75',
alpha = .5,
size = 1) +
scale_x_continuous(breaks = 1:9, limits = c(1,9)) +
xlab('prior draws') +
ggtitle('prior predictive check')## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning: Removed 165719 rows containing non-finite values (stat_bin).## Warning: Removed 2 rows containing missing values (geom_bar).
dataset1.brm.olr1 =
brm(dataset1.brm.olr.formula,
prior = dataset1.brm.olr.priors,
data = exp1.data,
backend = "cmdstanr",
file = "02-bayesian_stats/rds/dataset1.brm.olr1.rds"
)
summary(dataset1.brm.olr1)## Warning: There were 27 divergent transitions after warmup. Increasing
## adapt_delta above may help. See http://mc-stan.org/misc/warnings.html#divergent-
## transitions-after-warmup## Family: cumulative
## Links: mu = logit; disc = identity
## Formula: effectiveness ~ condition
## Data: exp1.data (Number of observations: 123)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept[1] -4.56 0.96 -6.78 -2.99 1.00 1070 513
## Intercept[2] -3.99 0.78 -5.82 -2.71 1.00 1583 1198
## Intercept[3] -3.29 0.56 -4.55 -2.32 1.00 2496 2468
## Intercept[4] -2.14 0.34 -2.85 -1.50 1.00 3348 3016
## Intercept[5] -1.26 0.26 -1.79 -0.76 1.00 3758 3330
## Intercept[6] -0.27 0.23 -0.74 0.18 1.00 3652 3366
## Intercept[7] 1.26 0.26 0.77 1.81 1.00 3532 3334
## Intercept[8] 2.22 0.33 1.62 2.91 1.00 3937 2992
## conditionno_graph 0.18 0.31 -0.40 0.81 1.00 3434 2673
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc 1.00 0.00 1.00 1.00 1.00 4000 4000
##
## Samples were drawn using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).dataset1.brm.olr2 =
brm(dataset1.brm.olr.formula,
prior = dataset1.brm.olr.priors,
data = exp1.data,
backend = "cmdstanr",
warmup = 1500,
iter = 2500,
control = list(adapt_delta = 0.99, max_treedepth = 15),
file = "02-bayesian_stats/rds/dataset1.brm.olr2.rds"
)
summary(dataset1.brm.olr2)## Family: cumulative
## Links: mu = logit; disc = identity
## Formula: effectiveness ~ condition
## Data: exp1.data (Number of observations: 123)
## Samples: 4 chains, each with iter = 2500; warmup = 1500; thin = 1;
## total post-warmup samples = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept[1] -4.58 1.00 -6.88 -3.00 1.00 2194 1948
## Intercept[2] -3.99 0.78 -5.76 -2.71 1.00 2806 2384
## Intercept[3] -3.30 0.55 -4.48 -2.32 1.00 3705 2806
## Intercept[4] -2.15 0.34 -2.85 -1.51 1.00 4378 3270
## Intercept[5] -1.26 0.26 -1.78 -0.76 1.00 4525 3093
## Intercept[6] -0.27 0.24 -0.73 0.19 1.00 3935 3132
## Intercept[7] 1.27 0.26 0.76 1.79 1.00 4089 3216
## Intercept[8] 2.23 0.34 1.58 2.92 1.00 4890 3122
## conditionno_graph 0.19 0.31 -0.43 0.82 1.00 4856 2853
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc 1.00 0.00 1.00 1.00 1.00 4000 4000
##
## Samples were drawn using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).plot(dataset1.brm.olr2)
dataset1.olr.y <- exp1.data$effectiveness
dataset1.olr.yrep <- posterior_predict(dataset1.brm.olr2)
ppc_hist(y = dataset1.olr.y,
yrep = dataset1.olr.yrep[1000:1007,], binwidth = .5)
ppc_dens_overlay(y = dataset1.olr.y,
yrep = dataset1.olr.yrep[2000:2030,])
ppc_stat_grouped(y = dataset1.olr.y,
yrep = dataset1.olr.yrep, group = exp1.data$condition)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
dataset1.olr.posterior_epred <-
tibble(condition = c('graph', 'no_graph')) %>%
add_epred_draws(dataset1.brm.olr2,
re_formula = NA,
allow_new_levels = FALSE) %>%
ungroup()
kable(head(dataset1.olr.posterior_epred))| condition | .row | .chain | .iteration | .draw | .category | .epred |
|---|---|---|---|---|---|---|
| graph | 1 | NA | NA | 1 | 1 | 0.0045580 |
| graph | 1 | NA | NA | 2 | 1 | 0.0063148 |
| graph | 1 | NA | NA | 3 | 1 | 0.0069149 |
| graph | 1 | NA | NA | 4 | 1 | 0.0391315 |
| graph | 1 | NA | NA | 5 | 1 | 0.0081217 |
| graph | 1 | NA | NA | 6 | 1 | 0.0248435 |
dataset1.olr.posterior_comparison <-
dataset1.olr.posterior_epred %>%
select(-c(.chain, .iteration, .row)) %>%
group_by(.category) %>%
compare_levels(variable = .epred, by = condition)
kable(head(dataset1.olr.posterior_comparison %>%
mean_qi()))| .category | condition | .epred | .lower | .upper | .width | .point | .interval |
|---|---|---|---|---|---|---|---|
| 1 | no_graph - graph | -0.0025112 | -0.0174983 | 0.0069861 | 0.95 | mean | qi |
| 2 | no_graph - graph | -0.0013803 | -0.0110381 | 0.0040049 | 0.95 | mean | qi |
| 3 | no_graph - graph | -0.0026421 | -0.0162269 | 0.0074813 | 0.95 | mean | qi |
| 4 | no_graph - graph | -0.0102726 | -0.0506452 | 0.0227708 | 0.95 | mean | qi |
| 5 | no_graph - graph | -0.0137966 | -0.0621373 | 0.0329472 | 0.95 | mean | qi |
| 6 | no_graph - graph | -0.0145299 | -0.0680404 | 0.0330627 | 0.95 | mean | qi |
dataset1.olr.posterior_comparison %>%
mean_qi(.epred) %>%
ggplot() +
geom_point(aes(y = .epred, x = .category), size = 3, color = DARK_PURPLE) +
geom_errorbar(
aes(ymin = .lower, ymax = .upper, x = .category),
width = 0,
size = 2,
color = DARK_PURPLE,
alpha = .5
) +
geom_hline(aes(yintercept = 0), linetype = "dashed", color = "gray") +
xlab('') + ylab('no_graph - graph') +
theme(axis.title.y = element_text(angle = 0, vjust = .5))
# Dataset 2
Load the data
dataset2 = readr::read_csv("02-bayesian_stats/data/exp2.csv") %>%
mutate(condition = condition == 'expansive') %>%
group_by(participant)##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## participant = col_double(),
## adjP10 = col_double(),
## adjP20 = col_double(),
## adjP30 = col_double(),
## adjPDiff = col_double(),
## condition = col_character(),
## gender = col_character(),
## index = col_double(),
## change = col_double(),
## subgroup = col_character(),
## orig = col_double(),
## BIS = col_character()
## )kable(head(dataset2))| participant | adjP10 | adjP20 | adjP30 | adjPDiff | condition | gender | index | change | subgroup | orig | BIS |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | 56.50000 | 55.80000 | 58.40000 | 1.90000 | TRUE | female | 69 | 3.362832 | expansiveHigh | 56.90000 | high |
| 2 | 27.87500 | 40.14286 | 36.00000 | 8.12500 | FALSE | female | 60 | 29.147982 | constrictiveLow | 34.67262 | low |
| 3 | 61.00000 | 59.00000 | 76.50000 | 15.50000 | TRUE | female | 55 | 25.409836 | expansiveLow | 65.50000 | low |
| 4 | 24.57143 | 32.75000 | 37.85714 | 13.28571 | FALSE | female | 79 | 54.069767 | constrictiveHigh | 31.72619 | high |
| 5 | 64.71429 | 54.00000 | 41.00000 | -23.71429 | TRUE | female | 69 | -36.644592 | expansiveHigh | 53.23810 | high |
| 6 | 35.14286 | 52.40000 | 45.60000 | 10.45714 | FALSE | female | 43 | 29.756098 | constrictiveLow | 44.38095 | low |
dataset2 %>%
mutate(c = as.factor(condition)) %>%
ggplot(aes(x = change)) +
geom_histogram(aes(y = ..density..),
binwidth = 6,
fill = DARK_PURPLE,
color = 'white') +
geom_density(size = 1, adjust = 1.5, color = 'lightgray') +
scale_x_continuous(limits = c(-200, 200)) ## Warning: Removed 1 rows containing non-finite values (stat_bin).## Warning: Removed 1 rows containing non-finite values (stat_density).## Warning: Removed 2 rows containing missing values (geom_bar).
Model 1
The first model is the BEST test model as described by Kruschke in the paper Bayesian estimation supersedes the t-test. In this model, \(\beta\) indicates the mean difference in the outcome variable between the two groups (in this case, the percent change in the BART scores). We fit different priors on \(\beta\) and set different weights on these priors to obtain our posterior estimate.
\[ \begin{align} y_{i} &\sim \mathrm{T}(\nu, \mu, \sigma) \\ \mu &= \alpha_{0} + \beta * x_i \\ \sigma &= \sigma_{a} + \sigma_{b}*x_i \\ \beta &\sim \mathrm{N}(\mu_{0}, \sigma_{0}) \\ \sigma_a, \sigma_b &\sim \mathrm{Cauchy}(0, 2) \\ \nu &\sim \mathrm{exp}(1/30) \end{align} \]
dataset2.brm.student.formula <- bf(
change ~ condition,
sigma ~ condition,
family = student()
)
kable(head(as_tibble(get_prior(dataset2.brm.student.formula, data = dataset2))))| prior | class | coef | group | resp | dpar | nlpar | bound | source |
|---|---|---|---|---|---|---|---|---|
| b | default | |||||||
| b | conditionTRUE | default | ||||||
| student_t(3, 22.6, 38.8) | Intercept | default | ||||||
| gamma(2, 0.1) | nu | default | ||||||
| b | sigma | default | ||||||
| b | conditionTRUE | sigma | default |
dataset2.brm.student.priorchecks = brm(
dataset2.brm.student.formula,
prior = c(
prior(normal(0, 2), class = "b"),
prior(cauchy(0, 2), class = "b", dpar = "sigma"),
prior(exponential(0.033), class = "nu"),
prior(student_t(3, 0, 5), class = "Intercept"),
prior(student_t(3, 0, 2), class = "Intercept", dpar = "sigma")
),
data = dataset2,
backend = "cmdstanr",
sample_prior = 'only',
file = "02-bayesian_stats/rds/dataset2.brm.student.priorchecks.rds"
)
# dataset2.student.yprior <-
# posterior_predict(dataset2.brm.student.priorchecks)
# ggplot() +
# geom_density(
# aes(x = dataset2.student.yprior),
# color = '#351c75',
# alpha = .5,
# size = 1
# ) +
# xlab('prior draws') +
# ggtitle('prior preditive check')dataset2.brm.student = brm(
bf(change ~ condition,
sigma ~ condition,
family = student()),
prior = c(
prior(normal(0, 2), class = "b"),
prior(cauchy(0, 2), class = "b", dpar = "sigma"),
prior(exponential(.033), class = "nu"),
prior(student_t(3, 0, 5), class = "Intercept"),
prior(student_t(3, 0, 2), class = "Intercept", dpar = "sigma")
),
data = dataset2,
backend = "cmdstanr",
file = "02-bayesian_stats/rds/dataset2.brm.student.rds"
)
dataset2.brm.student$prior## prior class coef group resp dpar nlpar bound
## normal(0, 2) b
## normal(0, 2) b conditionTRUE
## cauchy(0, 2) b sigma
## cauchy(0, 2) b conditionTRUE sigma
## student_t(3, 0, 5) Intercept
## student_t(3, 0, 2) Intercept sigma
## exponential(0.033) nu
## source
## user
## (vectorized)
## user
## (vectorized)
## user
## user
## usersummary(dataset2.brm.student)## Family: student
## Links: mu = identity; sigma = log; nu = identity
## Formula: change ~ condition
## sigma ~ condition
## Data: dataset2 (Number of observations: 80)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 21.15 4.95 11.41 31.04 1.00 2751 2071
## sigma_Intercept 3.37 0.20 2.98 3.76 1.00 2331 2293
## conditionTRUE -0.15 1.92 -4.01 3.46 1.00 3711 2961
## sigma_conditionTRUE 0.14 0.22 -0.31 0.55 1.00 3648 2440
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## nu 4.87 5.87 1.59 16.82 1.00 1781 1303
##
## Samples were drawn using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).plot(dataset2.brm.student)
dataset2.student.y <- dataset2.brm.student$data$change
dataset2.student.yrep <- posterior_predict(dataset2.brm.student)
dataset2.student.epred <- tibble(condition = c(TRUE, FALSE)) %>%
add_epred_draws(dataset2.brm.student, re_formula = NA, allow_new_levels = TRUE) %>%
ungroup()ppc_hist(y = dataset2.student.y,
yrep = dataset2.student.yrep[100:107,], binwidth = 10)
ppc_dens_overlay(y = dataset2.student.y,
yrep = dataset2.student.yrep[2000:2030,])
ppc_stat_grouped(y = dataset2.student.y,
yrep = dataset2.student.yrep,
group = dataset2$condition,
binwidth = 5)
dataset2.student.epred_comparison <- tibble(condition = c(TRUE, FALSE)) %>%
add_epred_draws(dataset2.brm.student, re_formula = NA, allow_new_levels = FALSE) %>%
ungroup() %>%
select(-c(.chain, .iteration, .row)) %>%
compare_levels(variable = .epred, by = condition) %>%
rename(diff = .epred)
kable(head(dataset2.student.epred_comparison))| .draw | condition | diff |
|---|---|---|
| 1 | TRUE - FALSE | -2.918310 |
| 2 | TRUE - FALSE | 2.367050 |
| 3 | TRUE - FALSE | 3.045840 |
| 4 | TRUE - FALSE | -0.455904 |
| 5 | TRUE - FALSE | -1.398820 |
| 6 | TRUE - FALSE | 2.193080 |
dataset2.student.epred_comparison %>%
select(diff) %>%
mean_qi() %>%
ggplot() +
geom_point(aes(x = diff, y = condition), size = 3, color = DARK_PURPLE) +
geom_errorbarh(
aes(xmin = .lower, xmax = .upper, y = condition),
height = 0,
color = DARK_PURPLE,
alpha = .5,
size = 2
) +
geom_vline(aes(xintercept = 0), linetype = "dashed", color = "gray") +
coord_cartesian(ylim = c(0, 2), xlim = c(-5, 5)) +
scale_y_discrete(label = c('expansive - not expansive')) +
xlab('') + ylab('')## Adding missing grouping variables: `condition`